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§ ; ABSTRACT 

' The r-modes of neutron stars can be driven unstable by gravitational radiation. While linear 

r T ^ perturbation theory predicts the existence of this instability, linear theory can't provide any 

information about the nonlinear development of the instability. The subject of this paper is the 
weakly nonlinear regime of fluid dynamics. In the weakly nonlinear regime, the nonlinear fluid 
equations are approximated by an infinite set of oscillators which are coupled together so that 
terms quadratic in the mode amplitudes are kept in the equations of motion. In this paper, 
the coupling coefficients between the r-modes are computed. The stellar model assumed is a 
polytropic model where a source of buoyancy is included so that the Schwarzschild discriminant 
is nonzero. The properties of these coupling coefficients and the types of resonances possible are 
discussed in this paper. It is shown that no exact resonance involving the unstable I = rn = 2 
r-mode occur and that only a small number of modes have a dimensionless coupling constant 
larger than unity. However, an infinite number of resonant mode triplets exist which couple 
indirectly to the unstable r-mode. All couplings in this paper involve the I > m r- modes which 
only exist if the star is slowly rotating. This work is complementary to that of Schenk et al. 
(2002) who consider rapidly rotating stars which are neutral to convection. 

Subject headings: instabilities — stars: neutron — stars: oscillations — stars: rotation 
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1. Introduction 

Gravitational radiation can drive nonaxisymmetric normal modes of rotating star unstable (Chan- 
drasekhar 1970; Friedman & Schutz 1978b; Friedman 1978) (CFS). The CFS instability affects nonaxisym- 
metric modes of stars which counter-rotate when viewed in a reference frame rotating with the star, but co- 
rotate when viewed in the inertial frame. The nonaxisymmetric modes are predicted to radiate gravitational 
radiation and the emission of radiation drives the instability. This provides a mechanism for converting the 
star's rotational energy into gravitational radiation, which will cause the star to spin down to slower angular 
velocities while providing a strong source of gravitational radiation which could potentially be detected. 

The traditional analysis of the CFS instability has been in terms of the instability of the fundamental 
pressure modes (f-modes) of rapidly rotating neutron stars (eg. Stergioulas 1998). However, it was suggested 
by Andersson (1998) that fluid modes of neutron stars driven by the Coriolis force (Rossby waves, or r- 
modes) can also be driven unstable by the CFS mechanism. It has since been shown that the r-modes are 
indeed unstable within the context of linear perturbation theory (Friedman & Morsink 1998) and that the 
growth times are short enough to be of astrophysical interest (Lindblom, Owen & Morsink 1998; Andersson, 
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Kokkotas & Schutz 1999). A review of recent results on the astrophysical relevance of the r-mode instability 
has been given by Andersson & Kokkotas (2001). 

The use of linear perturbation theory is sufficient to show that the modes are unstable, however a linear 
analysis can not describe the evolution of the instability when the amplitude of the perturbation grows 
large. It is inevitable, then, that linear perturbation theory can not be a sufficient description of an unstable 
system. In order to fully describe the growth and saturation of an instability, it is necessary to incorporate 
the nonlinear aspects of the system, either through a full numerical evolution of the nonlinear equations 
describing the system, or through a higher order perturbation expansion. 

One of the most interesting aspects of the r-mode instability to gravitational radiation reaction, is the 
prediction (Owen et. al. 1998) that the gravitational radiation could be detected with an advanced version 
of LIGO. However, the analysis made by Owen et. al. (1998) is dependent on the maximum amplitude 
the unstable mode is assumed to have before nonlinear effects dominate. The value of the maximum mode 
amplitude is inaccessible from a purely linear perturbation analysis and it is necessary to consider the 
nonlinearities in the fluid equations of motion. First steps towards including nonlinearities in the study 
of r-modes have been made. Two numerical evolutions of the fully nonlinear fluid equations have recently 
been made. The first evolution (Stergioulas & Font 2001) is fully relativistic, however, changes in the 
gravitational field induced by the perturbation had to be neglected in order to make the computation feasible. 
This evolution found that over the course of 25 r-mode pulsation periods, no saturation was observed. A 
second numerical evolution (Lindblom, Tohline and Vallisneri 2001a, 2001b) was performed using Newtonian 
gravity and post-Newtonian gravitational radiation reaction terms and did not show any signs of amplitude 
saturation due to nonlinearities. Neither numerical evolution can definitively explain the lack of nonlinear 
saturation. 

An alternative approach to the nonlinear evolution problem is to consider the weakly nonlinear regime 
and to do a higher order perturbation analysis. This approach involves taking the general Hamiltonian 
for a rotating fluid and perturbing around the equilibrium. This allows the explicit evaluation of coupling 
coefficients between different modes, so that it is possible to make predictions about the maximum amplitudes 
allowed before the perturbation analysis breaks down. This paper will apply this approach to the nonlinear 
coupling between the r-modes. This paper is a complementary study to that of Schenk et al. (2002) who 
have derived the Hamiltonian for nonlinear perturbations of rotating Newtonian stars. Schenk et al. (2002) 
have explicitly evaluated the lowest order nonlinear coupling coefficients for stars which have a small ratio of 
N/Cl, where N is the Brunt- Vaisala frequency (due to buoyancy) and is the star's angular velocity. In this 
paper, the opposite limt is considered. This is equivalent to a slow-rotation approximation, while the work 
of Schenk et al. (2002) corresponds to rapid rotation. The results of numerical evolutions of the interactions 
between the coupled modes considered in this paper and those considered by Schenk et al. (2002) will be 
presented in a following paper (Arras et. al. 2002). 

The weakly nonlinear regime has been investigated by other authors for other types of modes and stars. 
The role of parametric resonances in the nonlinear coupling of modes of nonrotating stars were discussed 
in detail by Dzicmbowski (1982). This work showed how stable equilibrium solutions can lead to mode 
saturation and was extended to describe the behaviour of Cepheids by Dziembowski & Kovacs (1984). The 
role of nonlinear interactions between p-modes in the Sun were investigated by Kumar & Goldreich (1989) and 
Kumar, Goldreich, & Kerswell (1994). Amplitude saturation of g-modes in white dwarfs through nonlinear 
interactions was shown in a paper by Wu & Goldreich (2001). All of the above calculations have involved 
only the lowest order nonlinear terms. Van Hoolst (1994) has constructed a higher order Hamiltonian which 
includes terms which are quartic in the displacement. Higher order terms have also been included in a 
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relativistic study of radial oscillations (Sperhake, Papadopoulos, & Andersson 2001). None of the above 
studies has considered rotating stars. However, in order to be able to study the saturation of an unstable 
r-mode, it is necessary to include the complicating effect of rotation. In geophysics, r-modes in oceans are 
of great interest and nonlinear interactions between the r-modes in oceans have been studied by Longuet- 
Higgins & Gill (1967), Ripa (1981), and Pokhotclov ct al. (1995). However, these calculations are typically 
done in the "/3-plane" approximation where locally the surface of a sphere is replaced by a plane. This 
approximation does not describe very well the situation of interest in stellar physics, where the star as a 
whole is expected to oscillate. In the calculations where the curvature of the ocean is not neglected, the 
radial dependence of the mode's cigenfunction is ignored, and this is not valid for stellar physics. As already 
discussed above, Schcnk ct al. (2002) are the first to study the effect of nonlinear interactions amongst the 
modes of rotating stars. 

The structure of this paper is as follows. In section 2 linear perturbation theory is briefly reviewed. In 
section 3 an alternative form of the third order Hamiltonian is derived. In section 4 the form of the nonlinear 
equations of motion for oscillations of a rotating star are reviewed. The type of resonances allowed in the 
coupled r-mode system are discussed in section 5. Numerical results are presented in section 6 and their 
implications are discussed in the conclusion. 



The equilibrium state of a rotating star is described by its pressure p, density p, velocity v a and 
gravitational field In this work, only rigidly rotating stars are considered, so that the star's velocity field 
is of the form v a = Q(f> a where f2 is a constant angular velocity. Small perturbations about the equilibrium 
star are described by the displacement vector field £ a which connects a fluid element in the equilibrium star 
to a fluid element in the perturbed star. 

The Hamiltonian for the perturbation theory of a nonrclativistic fluid is a series expansion in the fluid 
displacement field of the form 



where each term in the series depends on the equilibrium fluid variables and £. The subscript on each term 
denotes the order at which the displacement field £ a appears in the term. If all terms up to and including 
the term 5H n are kept in the expansion, then the resulting equations of motion will be of order (n — 1) in 
£ a . In this notation, retaining only the second order Hamiltonian term results in linear perturbation theory. 

The second order Hamiltonian for a nonrelativistic rigidly rotating fluid is (Friedman & Schutz 1978a) 



2. Linear Perturbation Theory of Rigidly Rotating Nonrelativistic Fluids 



sH = 6H 2 (e,n+sH 3 (e,e,e)+ - +«/„(r,...,n + o(i£r +i ), 



(i) 




(2) 



where the Eulerian perturbations of the density and pressure are given by 



Sp 
Sp 



(3) 
(4) 



and the perturbed gravitational potential satisfies 



V 2 <5$ = AnGSp. 



(5) 
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The expression (2) can be physically interpreted as the energy of the perturbation as measured in the rotating 
frame and has also been named the rotating frame canonical energy by Friedman & Schutz (1978a), where 
its properties have been discussed extensively. 



3. An Alternative Expression for the Third Order Hamiltonian 

The expansion of the Hamiltonian up to and including terms third order in the displacement £ has been 
found by Kumar & Goldreich (1989) for the case of static fluids 1 . When the equilibrium star is rotating, the 
contributions to the perturbed Hamiltonian through the perturbation's kinetic energy is altered. Surprisingly, 
Schenk et al. (2002) have shown that the contribution from the kinetic energy of the perturbation is exactly 
second order in the displacement. As a result, all terms in the perturbation expansion of the Hamiltonian 
for a rotating star which are third order in the displacement (or higher order) are exactly the same as the 
terms in the perturbation expansion for a nonrotating star. This is an important result and simplifies the 
study of perturbations of rotating stars greatly. The third order term in the perturbation expansion of the 
Hamiltonian is 

SH 3 = 1 J dV [-2 p V a £"V b rV c C + p rC 6 rV a V fe V c $ (6) 

-3 p (Tx - 1) V • £ V a i; b V b e - p (V i - I) 2 + (V • o 3 

+3 P ee'V J V J M>] , 

as derived by Kumar & Goldreich (1989) and corrected by Schenk et al. (2002). 

Although the third order Hamiltonian (6) describes the nonlinear interactions of modes, it is not simple 
to see from this expression what the relative size of the interaction term is to other physical energy scales. 
For this reason, an alternative form of equation (6) will now be presented which is useful for finding the size 
of the nonlinear coupling between different r-modes. 

The goal is to replace the first term in (6) with a simpler expression by integrating by parts. The explicit 
calculation can be found in the Appendix. The final result is 

-2 J dVp V £ 6 V 6 £ c V c r = J dV [reTVaV fc V c p + 3(V • £) £ 6 £ c V c V 6 p 

-3V • (K)V a £ b V fc £ Q + 2 (V • 2 Z ■ Vp 

-2p (V • TV a (V • 0] • (7) 

While the expression (7) doesn't appear to be a simplification, it is useful, since the first term involves three 
derivatives of the equilibrium pressure. Recall that the third order Hamiltonian (6) has a term involving 
three derivatives of the equilibrium gravitational potential. This suggests that these terms be combined 
using the equilibrium equations of motion. The equation of hydrostatic equilibrium for a rotating star is 

iv a (« 2 ) + iv a p + V a $ = 0. (8) 
z p 



1 The terms involving the perturbed gravitational field given by Kumar & Goldreich (1989) arc incorrect and have been 
corrected by Schenk et al. (2002). 
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The strategy will be to take two covariant derivatives of the equation of hydrostatic equilibrium and use the 
result to simplify the third order Hamiltonian. 

In the case of rigid rotation, 

V a V b V c (« 2 ) = 0. (9) 
To prove equation (9), use Cartesian coordinates, so that v 2 = Q 2 (x 2 + y 2 ) and V a = d a . Then 

V a {v 2 ) = 2n 2 (xd a x + yd a y) (10) 
V 6 V a (« 2 ) = 2n 2 (d a xd b x + d a yd b y) (11) 

from which eq. (9) follows. 

Taking the two derivatives of the equation of hydrostatic equilibrium and taking the inner product with 
three displacements, we find the following equality 

-rCT ( ^^-VapVbpVcP + -^V a pV b pV cl - — V a \7 bP V c p) . (12) 

V pi p pi p ipp ) 

The hydrostatic equilibrium equation (12) and the integration by parts formula (7) can now be substi- 
tuted into the equation for the third order Hamiltonian (6) resulting in the equivalent third order Hamiltonian 



sh 3 = i J dv s P (v a ev b e - ±-^ev a v bP ^ 



(^ T ^-VapV bP V c p+ -^V aP V &P V c7 

V pi v pi p 

+2(v-o ((v-o^-vp-pe-v(v-o) 

-3 p (r x - 1) V • % V a £ b V h C - p ((T, - l) 2 + ) (V • o 3 



d In p, 

+3 p ee'ViVj(5$] . (13) 

Equation (13) was derived using only the assumption of rigid rotation. No assumption about the type of 
perturbation has been used. 

It is now of interest to consider the coupling of r-modes at third order in perturbation theory. The 
traditional r-modes are solutions of the linearized Euler equation in the limit of small rotation rate, Q = 
ily^R 3 /M < 1 and the slow rotation approximation will be used for the remainder of this paper. In the slow 
rotation approximation, the r-mode displacement satisfies (eg. Friedman & Morsink 1998) 

V-f~5p~£- Vp~£- V 7 -(9(n 2 ), (14) 

while in Appendix B it is shown that the terms 

v a e^ b c ~ rc b v a v bP ~ o(i) (15) 

in the slow rotation expansion. 

In the slow rotation approximation, adopting the Cowling approximation the third order Hamiltonian 
coupling r-modes together is 



5Hi 



= ±JdVS P (v a ev b e - i-^^V„V 6 p) + 0(ti 4 ). (16) 
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Note that the third order Hamiltonian is 0(Vl) 2 in the slow rotation expansion. (The term dropped by 
assuming the Cowling approximation is also 0(&) 2 .) This has important consequences, since the energy of 
the r-modes (the term SH 2 in the perturbation expansion) is also 0{Cl) 2 (Friedman & Morsink 1998). If 
the third order term was first order in the dimensionless angular velocity, then the third order term would 
dominate over the second order terms whenever the amplitude of the r-mode grew larger than the square of 
the angular velocity. Instead, we see that the third order term in the Hamiltonian is of the same order (in the 
small angular velocity limit) as the second order term. As a result, if the third order terms are responsible 
for saturation, the limit on the r-mode amplitude which they set will be independent of angular velocity. 
However, this is strictly only true if no dissipation is present. Since the dissipative timescales will depend 
on angular velocity, there will always be an angular velocity dependence in the saturation amplitude. 

In order to investigate the possible saturation of an unstable r-mode through couplings with other r- 
modes at third order, it is necessary to solve for the r-mode at a high enough order in the slow rotation 
approximation that the pressure perturbation is non-zero. For if only the lowest order in angular velocity 
terms are kept, the r-mode pressure perturbation will vanish and the third order Hamiltonian will also vanish 
(Schenk et al. 2002). 



4. Nonlinear Equations of Motion 

The equations of motion for a general perturbation at any order of perturbation theory can be found 
in a straight-forward manner if the equilibrium star is nonrotating. The general procedure is to expand the 
displacement into a sum over normal mode solutions 

^ = E c 4*)(aW (17) 

A 

where Q A (x)e~ lulAt is a solution of the linearised equations of motion. The expansion coefficients, c A (t) 
are found by substituting in the expansion (17) into Euler's equation and making use of the orthogonality 
properties of the functions Q A . This procedure is not as straight- forward when the background star is 
rotating. The main complication is that the normal mode solutions are not orthogonal with respect to the 
usual inner product, 

< Ca, (b >= J dV pC a Cb ^0 for A^B. (18) 

If such a procedure were to be followed, the resulting equations of motion would be coupled at linear order 
(Schenk et al. 2002). The coupling of equations at linear order can be circumvented by using a phase space 
expansion of the displacement, which combines equation (17) with the expansion (Schenk et al. 2002) 

i = J2(-i"A)c A (t)(A(x). (19) 

A 

The resulting equations of motion derived by Schenk et al. (2002) for the expansion coefficients (including 
third order terms in the Hamiltonian) are 

c A (t) + iu A c A {t) = ^Y1 — c B(t)ch(t) (20) 

where e A is the mode's energy in the rotating frame at unit amplitude defined by 

5H 2 = ^e A |c4| 2 (21) 
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and the nonlinear coupling coefficient k is denned by 

5H 3 = -\j2 CA{t)c B {t)c c {t)K ABC . (22) 

ABC 

There is freedom to choose the amplitudes of the spatial mode functions £• One simple scheme is to 
choose these amplitudes so that the rotating frame energies all take the value 

e A = MR 2 n 2 (23) 

which is the same order of magnitude as the star's kinetic energy. With this choice, \ca\ 2 is the ratio of the 
mode's energy to the star's energy. Since the coupling coefficients kabc have units of energy, the fraction 
kabc / e A corresponds to a dimensionless fraction of the star's energy. When the mode amplitudes satisfy 

\cb\\cc\ _ e A ^ 
\ca\ kabc 

the nonlinear terms in the equation of motion (20) will dominate over the linear terms, which signals the 
breakdown of the weakly nonlinear regime. If (24) is satisfied, it will be necessary to include higher order 
nonlinear terms in the equations of motion. 



4.1. Coupling between r- modes 

In the case of the r-modes, the coupling coefficients are 

kabc = J dV 



Spa [V a ( B V b (g,-—( B C b cVaV bP \ (25) 



PI ' 

+5 PB (v a ( b c V b ( a A - — ( A ( c V a V bP ) +S PC (v„CaV 6 Cb - — aCBV a V 6 p 
V pi / V pi 

where Spa = — I\pV • (a — (a • Vp. The coupling coefficients for the case of three r-modes are of second order 
in angular velocity. In the case of r-modc couplings, a dimensionless coupling coefficient, kabc is defined by 

kabc ,„m 
k A bc = j^^, (26) 

so that the magnitude of the kabc coefficients are independent of the star's angular velocity. This choice of 
dimensionless coupling coefficient also has the nice feature that the coefficients are independent of the star's 
mass, once a polytropic model has been chosen. 

The r-mode frequencies, given quantum numbers I a and m^, are 

oja = , 2 , r ! lAn ^ +C(lA,mA,k A )nn 2 , (27) 

where the frequency correction C(Ia, tua, ^a) depends on Ua, the number of radial nodes in the cigenfunction, 
as well as the equation of state and buoyancy law. It is convenient to introduce dimensionless frequencies 
Co A are defined by 

* = (28) 
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With the new dimensionless coefficients and frequencies, the equations of motion coupling r-modes 
together are 

c A {t) + iu A c A {t) = lAC A (t) + V K* ABC c* B (t)c* c (t) (29) 

1 BC 

where the time coordinate is now measured in units of the star's spin period, and an external damping/driving 
term proportional to "/ A has been included. When the external damping and driving terms are neglected, 
equation (29) is independent of angular velocity. As a result, as the star spins down, the nonlinear coupling 
equations are unchanged if no dissipation is present. When dissipation is included, the strength of the 
dissipation terms changes with time if the star spins down. The scalings with angular velocity which have 
been used are only strictly correct in the slow rotation approximation. However, calculations of r-modes 
for rapidly rotating stars (Lindblom & Ipser 1999; Yoshida et al. 2000; Karino et al. 2000) have shown that 
these scalings are not a bad approximation when the star is rapidly rotating. 



5. Resonance Condition and Selection Rules 

The coupling coefficients k A bc involve integrals of three r-mode eigenfunctions over all space. This 
leads to two selection rules on the possible values of I and m for each mode. Integration over (f> gives the 
selection rule on the azimuthal quantum numbers, 

m A + m B + m c = 0. (30) 

Since the integration is over all space, the overall parity of the integrand must be even. Since r-modes have 
axial parity, the parity of an r-mode with angular momentum quantum number I is (— It follows then 
that the only allowed couplings of r-modes must obey 

U + Ib + Ic + 1 = (mod 2). (31) 

It is simple to see that as a result of these selection rules, there are no couplings between r-modes of 
barotropic stars (ie. neutrally convective stars) in third order perturbation theory. Stars which are barotropic 
only have r-mode solutions with I = \m\. The I selection rule (31) in this case is 

\m A \ + \m B \ + \m c \ + 1 = (mod 2). (32) 

However it is a property of the integers that 

\m A \ + |ms| + \mc\ = m A + m B + m c (mod 2). (33) 

In order to have nonzero coupling coefficients between r-modes of barotropic stars, we would require 0=1 
(mod 2), which is clearly impossible (Schcnk et al. 2002). 

One might expect that the usual triangle inequality which occurs in the addition of angular momentum 
should hold here as well. However, since the formula for the coupling coefficients involves the pressure 
perturbation, the triangle inequality is modified. If the lowest-order axial term in the r-mode expansion 
corresponds to quantum number lc, the pressure perturbation will have quantum number lc± 1 (Saio 1982). 
The modified triangle inequality is then 



Ic-1<Ia + Ib and \l B - l A \ < l c + 1. 



(34) 
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In order to find nonzero coupling coefficients between r-modes, the perturbations must have a nonzero 
Schwarzschild discriminant so that r-modes with I ^ \m\ may exist. Including a Schwarzschild discriminant 
in the perturbation theory is equivalent to introducing a source of buoyancy. Some physical mechanisms for 
the inclusion of buoyancy in neutron stars which have been suggested are finite temperature (McDermott, 
Van Horn, & Scholl 1983), composition gradients (Rciscncggcr & Goldrcich 1992), and accretion (Bildsten & 
Cutler 1995). The general effect of buoyancy can be modelled by introducing an adiabatic index Ti different 
from the polytropic index 7 = 1 + 1/N. The dimcnsionless parameter 1 —-f/Ti is typically quite small for all 
of these buoyancy models. For instance, for the composition gradient buoyancy introduced by Reisenegger 
& Goldreich (1992), 

1 - -j- ~ 3 x 10~ 3 — (35) 

1 1 Pnuc 

in the core of the star, where p nuc — 2.8 x 10 14 gcm~ 3 . The computation of r-modes including a realistic 
adiabatic index which varies with position is a challenging problem, since the equations for the / ^ \m\ r- 
modes are not valid at any places where 1 — 7/1^ vanishes. At present, no calculation of r-modes of neutron 
stars with varying Ti exists in the literature. A reasonable compromise is a simpler model in which T 1 is a 
constant. The r-modes of polytropic stars with constant Ti have been examined in great detail by Yoshida 
& Lee (2000). Following the work of Yoshida & Lee (2000), we will also adopt similar buoyancy models 
in the present paper. An important property of the I ^ \m\ r-modes has been pointed out in the work of 
Yoshida & Lee (2000): the I ^ \m\ r-modes only exist in the limits of slow rotation compared to the break-up 
velocity (ft <C 1) and slow rotation compared to the Brunt- Vaisala frequency which places an approximate 
upper limit of Cl ~ yT — 7/T1 on the angular velocity. 

In order to make fully nonlinear numerical simulations computationally possible, the physics is generally 
simplified to the level of a perfect, zero-temperature fluid, as in the simulations by Stergioulas & Font (2001) 
and Lindblom et al. (2001a, b). For this reason, these simulations can not include any source of buoyancy 
and the stars can only support the oscillations of r-modes with I = \m\. Since the arguments presented above 
show that there are no nonlinear couplings between I = \m\ r-modes at third order in perturbation theory, 
the fully nonlinear codes do not access the same physical coupling mechanisms explored in the present paper. 
At linear level, the codes written by Stergioulas & Font (2001) and Lindblom et al. (2001a, b) model the 
perturbations known as hybrid modes (Lockitch & Friedman 1999). Nonlinear couplings between the hybrid 
modes have been computed by Schenk et al. (2002). The hybrid modes are good description of a star's 
low-frequency modes when the Brunt- Vaisala frequency is very small compared to the spin frequency. 

It is well known from the theory of nonlinear oscillations that whenever a resonant match of frequencies 
occurs in a system consisting of three coupled oscillators, it is possible to transfer energy from a large 
amplitude mode to a small amplitude mode. The condition for resonance is that the detuning, defined by 

Aluabc = u A + lu b + ujc, (36) 

should be close to zero. If any resonances occur amongst the direct couplings with the I = m = 2 r-mode, 
the unstable mode's amplitude could be limited by the coupling, even if the coupling coefficient is small. 

In the numerical simulations of Stergioulas & Font (2001) and Lindblom et al. (2001a, b), the initial 
amplitude of the unstable I = m = 2 r-modc is set to a small value and all other mode amplitudes are set 
to zero. In terms of the equations of weakly nonlinear system, this corresponds to the I = m = 2 r-modc 
acting a a "parent" or source for a later generation of modes. As an example, the parent mode will carry 
the label A = 0. Since all other modes initially have zero amplitude, the first generation of daughter modes 
will be only those with nonzero coupling coefficients «ioo, where the label 1 refers to an excited daughter 
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mode. In the case of the CFS instability of r- modes, the unstable parent mode is the I = \m\ = 2 r-mode. 
The daughter modes which can be excited through the third order coupling can be found using the selection 
rules (30) and (31), choosing l A = l B = 2, m,A = rriB = 2 and solving for lc and mc- The daughter modes 
which can be excited are shown in Table 1. The detuning can be easily calculated by using the formula 
for the r-modc frequency in the slow rotation limit f2 <C 1. For this first generation of couplings, with 
I a = m,A = Ib = m B = 2, and mc = —4, the detuning is 

Am abc = 4 _ 8 = 4(Z c + 3)(/ c -2) 

fi 3 l c {lc + l) 3l c {lc + l) ' 1 ' 

which is non-zero, since the selection rules require lc = 5. 

Similarly, a second generation of daughter modes can be excited through couplings of first generation 
daughters with the unstable r-mode or through the coupling of two daughter modes. The second generation 
daughter modes are also displayed in Table 1. It should be clear that modes excited at any generation must 
have an even value of m. Furthermore, it is impossible to generate certain modes, such as an r-mode with 
m = 2 and an odd value of I if the only nonzero initial amplitude corresponds to the unstable I = m = 2 
r-mode. This type of coupling, where the parent mode acts as a source for a daughter mode will be denoted 
a direct coupling. In the case of the second generation direct couplings, the detuning is either 

^ABC 2(Z 2 + i + 10) 

51(1 + 1) > (38) 

in the case of couplings involving Ib = 5 and = -4 or 

Aujabc _ 2 (7/ 2 + 7/ - 6) 



n i5i(z + i) 



(39) 



which is also nonzero for whole number values of I. This shows that a numerical simulation which sets to zero 
the initial values of all modes except for the I = m = 2 r-mode will only involve non-resonant interactions. 

In a realistic star it is expected that the other modes will have small but nonzero amplitudes. When 
is it possible to have a resonance in any triplet of r-modes including the I a = tha = 2 r-modc? Suppose 
we choose values lc — I and mc = m for the quantum numbers for mode C. From the m-selection rule 
nig = — (to + 2) and if we choose Ib > I the 1-selection rule and the triangle inequality restrict l B to the 
form l B = I + 2fc + 1 where k is either or 1. The general formula for the detuning in this case is 

Am abc = 2 _ 4 4m(2fc + l)(/ + fc + l) 

ft ~3 l B {l B + l) 1(1 + 1)I b {Ib + 1) ' ( ' 

In Figure 1 it can be seen that the detuning can't vanish if we are restricted to nonzero values of m and 
mB- Hence there are no direct resonant couplings with the I = to = 2 r-mode. However, the values of the 
detuning are typically smaller than unity and quickly converge to 2/3 for large values of I. 

Any three r-modes with an identical odd value of I will be resonant if they satisfy the m-selection rule. 
This occurs because the detuning for modes with I a = Ib = lc = I is 

^^ = T(iTT) {mA + mB + mc) (41) 

which vanishes when equation (30) is satisfied. The Z-selection rule (31) demands that I must be odd, so that 
these resonant triplets do no involve the I = m = 2 mode. Since the resonance criteria occurs for all odd 
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values of I, there are an infinite number of resonances, none of which couple directly with the I = m = 2 r- 
mode. These couplings are indirect resonant couplings. These resonances will not be exact, since the r-mode 
frequencies have the small correction terms shown in equation (27). In the slow rotation limit, however, 
these corrections are very small, and these triads will be very close to resonance. 

6. Numerical Results 

In this paper, the r-mode eigenfunctions of polytropic stars are solved by a numerical method, similar 
to the method described by Saio (1982), which keeps terms in the perturbed velocity that are third order 
in the angular velocity. The code's accuracy was checked against the results published by Saio (1982) for 
the case of modes where 7 ^ Ti and against the results of Lindblom, Mendell & Owen (1999) for the case 
of modes of stars with 7 = Ti. Good agreement was found in both cases. Once the r-mode eigenfunctions 
for large number of modes has been found, the coupling coefficients can be computed by directly integrating 
equation (25). In these calculations, all r- modes with I < 10 and less than 3 radial nodes were computed. 

There are two types of three- mode couplings which occur amongst the r- modes of buoyant stars. Triplets 
of modes which include the I — m = 2 r-mode are always nonresonant. These direct nonresonant mode 
couplings will be examined in section 6.1. Resonant mode couplings can occur triplets which do not include 
the I = m = 2 r-mode. This type of coupling is an indirect resonant coupling and will be examined in section 
6.2. 



6.1. Direct Nonresonant Mode Couplings 

Table 1 shows the dimcnsionlcss coupling coefficients for generating the first and second generations 
of r-modes. The equilibrium stellar model is a N — 1 polytrope and the adiabatic index Ti = 1.9. The 
N = 1 polytrope was chosen to allow the closest comparison with the nonlinear simulations of Stergioulas 
& Font (2001); Lindblom et al. (2001a, b) who all make use of the same polytropic model. However, since 
Ti 7^ 1 + 1/N, the types of couplings discussed in this paper will not be seen in the nonlinear simulations. 
The size of the coupling coefficients in Table 1 are all small. Equations (24) and (26) imply that the 
unstable r-mode would have to grow to an amplitude so that the fraction of energy in the mode to the star's 
energy is of order 1/k 2 before higher order nonlinear terms would have to be considered. Since the coupling 
coefficients are of order 10~ 2 , this suggests that the unstable r-mode would need to have an energy which is 
approximately 10 4 times the energy of the star before higher order nonlinear terms become important. In 
other words, the weakly nonlinear approximation holds for all physical values of mode energy. 

The dependence of the coupling coefficients on the equation of state and buoyancy law was also examined. 
The coupling coefficients corresponding to the first generation mode coupling (for modes with zero radial 
nodes) are shown in Table 2. There appears to be very little dependence on the equation of state or on 
buoyancy. The lack of dependence on buoyancy can be understood by considering equation (25) which defines 
the coupling coefficients. The coupling coefficients only depend on equilibrium quantities (the pressure) and 
the eigenfunctions of the perturbations. There is no explicit dependence in (25) on the eigenvalue C(l, m, k) 
defined in equation (27). As shown by Yoshida & Lee (2000), the eigenvalues and eigenfunctions of the I = \m\ 
r-modes have very little dependence on the buoyancy law. The only dependence on buoyancy can come 
through the I > \m\ r-modes. For the eigenvalues of these modes there is a strong dependence on buoyancy, 
which we have found varies approximately as C(l, m, k) ~ (Ti — 7) _1 . (This scaling can be also be deduced 
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from the form of the equations presented in the work of Saio (1982).) However, as long as the eigenvalue 
C(l, m, k) and the dimensionless angular velocity Cl are small enough that C(l, m, k)Vl 2 « 2m/l(l + 1), the 
form of the eigenfunctions are almost independent of the buoyancy. In order to satisfy this inequality for 
angular velocities as large as £l <~ 0.1, we require that Ti — 7 > 10~ 3 . As the difference Ti — 7 decreases to 
zero, the I > \m\ r-modes cease to exist (Yoshida & Lee 2000) so we would expect that at small values of 
Ti — 7 the character of these r-modes will be greatly modified and the values of the coupling coefficients will 
be altered. In Table 2 the values of the coupling coefficients for 7 = 2 are practically the same for all values 
of T\. It is only once T\ — 7 falls below 8 x 10~ 4 that any variation in the coupling coeffients is found. 

The largest coupling coefficients for modes which couple to the I a = rriA = 2 r-mode are shown in Table 
3. The largest coefficients appear to be those which couple (Zs,ms) = (I, —1), with (lc,mc) = (l + l, — 1). 
As long as these modes have a small but nonzero amplitude, it will be possible for these modes to be excited. 
In this case, these modes could be the most important type of interaction which could conceivably limit the 
amplitude of the unstable r-mode. 

As with all other types of modes of nonbarotropic stars, there are an infinite number of r-mode solutions 
for each value of I and m. Each solution can be labelled by the quantum number k > which corresponds 
to the number of radial nodes in the eigenfunction. The CFS instability in the r-modes is strongest in the 
fundamental k = modes, so only the fundamental / = m = 2 r-mode will be considered in this work. 
However, it is conceivable that the unstable mode could excite daughter modes with k > 0. In Table 
1, the dependence on number of nodes for the first generation daughter mode is shown. It is intriguing 
that the magnitude of the coefficients does not seem to fall off rapidly with increasing node number. This 
behaviour is due to the weighting of the integrals in the coupling coefficients, which favours the regions 
of the eigenfunctions closest to the surface of the star. Since most of the nodes in a k > eigenfunction 
occur in the middle regions of the star, and the eigenfunctions are peaked near the surface, there is not a 
large difference in the contributions of different /c-eigenfunctions to the coupling coefficient integrals. All of 
these higher-order eigenfunctions are able to participate in interactions with the I = m = 2 r-mode. This 
increases the number of modes which can act as a drain of the unstable mode's energy. In order to evaluate 
an upper bound on N, it is necessary to compute high-radial order r-modes, which becomes difficult since 
these modes oscillate rapidly. A more suitable approach would be an approximation scheme such as the 
WJKB method. This type of approach would allow an asymptotic calculation of the coupling coefficients 
for high-radial order modes. However, it seems unlikely that the high-radial order modes would participate 
in the type of interactions described above, for they would be rapidly damped by viscosity. In this case of 
a series of modes with approximately equal frequencies and coupling coefficients, the most slowly damped 
mode is the one which places the strongest limit on the unstable parent mode (Wu & Goldreich 2001). In 
this case, the fundamental mode is damped the least by viscosity, so it should be the most important. 



6.2. Indirect Resonant Couplings 

In section 5 it was shown that there are an infinite number of resonant triads which do not involve the 
I = m = 2 r-mode. These resonant triads only couple indirectly to the unstable r-mode. 

The coupling coefficients for resonant mode interactions are shown in Table 4. We note that the coupling 
coefficients in Table 4 are at least an order of magnitude larger than the coefficients in Table 1. Due to the 
infinite number of possible couplings of this type, only a small selection of triplets are shown in Table 4. 

In order to examine the possible effect of indirect resonant couplings on the amplitude of the I = m = 2 
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r-mode, we can consider a toy model consisting of only three modes. Mode A has quantum numbers 
I a = tua = 2. Mode B has quantum numbers Ib — 5 and tub — —4. Modes A and B can interact through 
a triplet involving two of mode A and one of mode B, with coupling coefficient fc^ := haab- The AAB 
triplet is a nonresonant interaction, for the reasons discussed in the previous subsection. The third mode, 
mode C, has quantum numbers lc — 5 and lc — 2. Mode C can interact with mode B through a triplet 
involving two of mode C and one of mode B, with coupling coefficient kc ■— kccb- The CCB coupling is 
resonant. Interactions between A and C are not allowed by the selection rules discussed at the beginning of 
this section. 

Since the CCB coupling is resonant, it is possible for an efficient transfer of energy to occur which would 
limit the amplitude of mode B. Suppose that mode B is the most important mode for saturating mode A 
when only triplets interacting with A are considered. As soon as mode C is included, the rapid transfer of 
energy between B and C could drastically alter the saturation of mode A. This type of situation is generic, 
since all allowed triplets involving the I = m = 2 r-mode also involve one r-mode with an odd value of I. 
Since r- modes with odd I always have resonances (except for / = 1), this type of indirect resonant interaction 
is possible. A similar situation was encountered in the study of g-modes by Wu & Goldreich (2001). In their 
study, Wu & Goldreich (2001) found that such couplings slowed down the transfer of energy between modes 
A and B, but were unable to halt the energy transfer. A similar study for r-modes will be the topic of a 
future investigation. 

7. Conclusions 

In this paper, we derived an alternative form of the Hamiltonian describing the interaction of r-modcs 
at lowest nonlinear order. This Hamiltonian was used to compute coupling coefficients between r-modes 
of polytropic stars with a source of buoyancy. In the case that only the I = m = 2 CFS unstable r-mode 
has a nonzero amplitude, none of the couplings are resonant, and the coupling coefficients are all small. If 
small but nonzero amplitudes are allowed for the damped r-modes, larger couplings can occur. These triads 
of modes with the largest coupling coefficients are the ones most likely to contribute to the saturation of 
the unstable r-mode's amplitude. An unusual aspect of this problem which we have pointed out is that 
there are an infinite number of resonant triads which couple indirectly to the unstable r-mode. It should be 
noted that all of the couplings considered in this paper involve the I > \m\ r-modes which only exist in the 
limit of slow rotation. This study is complementary to the work of Schenk et al. (2002) who examined the 
coupling coefficients for the modes of rotating stars which are neutral to convection, which is a very good 
approximation when the star is rapidly rotating. The results of a numerical evolution of the coupled system 
of oscillators considered in this paper and in the paper by Schenk et al. (2002) will be presented in a future 
paper (Arras et al. 2002). 

This paper has been concerned only with the couplings between different r-modes. However, it may 
be possible to couple r-modes with other modes, such as the g-modes or p-modes. The g-modes are highly 
distorted by rotation (Bildsten, Ushomirsky, & Cutler 1996), so the effect of their coupling with r-modes is 
difficult to compute in the slow-rotation approximation. Instead, for rapid rotation, the inertial or hybrid 
modes discussed by Lockitch & Friedman (1999), and Yoshida & Lee (2000) are a better approximation to 
the g-modes and r-modes. Hence the coupling between hybrid modes studied by Schenk et al. (2002) should 
be sufficient to study all of the different types of couplings between different families of modes. 

When the driving and damping of modes is included in the evolution equations describing the coupled 
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oscillators, it is possible for equilibrium solutions with constant amplitude to result (Dziembowski 1982; 
Dziembowski & Kovacs 1984; Wu & Goldreich 2001). This leads to a limiting amplitude on unstable 
modes. This limiting amplitude depends on the coupling coefficients, frequency detuning and the driving 
and damping rates of the modes. With the state of understanding of these parameters coming from this 
paper and the work of Schenk et al. (2002), it will be possible for numerical simulations of mode-mode 
coupling to show under what conditions it is possible for saturation of the unstable r-mode to occur. 
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research was supported by the Natural Sciences and Engineering Research Council of Canada and by NSF 
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appears. The purpose of this appendix is to transform the integral (Al) into a form more suitable for 
evaluating the r-mode coupling coefficients. The goal is to transfer the derivatives acting on the displacements 
to derivatives acting on the pressure. This can be done by integrating by parts a many times. All surface 
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A. Integration by Parts 



In equation (6) the integral 




(Al) 



- 16 - 



terms are proportional to either the pressure or the radial derivative of the pressure, both of which vanish 
on the surface of the star. 

Integrating by parts, and throwing away the surface terms, 

I = - f dVCV c ( P V a eV b e) =Ia + Ib + Ic, (A2) 

where 

I a = - J dVV c pCV a e V b f (A3) 
Ib = - J dVpCV a V c eV b C (A4) 



Ic 



JdV P CV a eV b (V-0- (A5) 



The integral I a has one derivative of pressure appearing. Integrating by parts once more will produce 
a term involving two derivatives of pressure, since 

Ia = f dVZ c V b (V cP CV a e) = Iai + Ia2 + Ia3, (A6) 

where 

Iai = J dVli c V c V b p C V a £ 6 (A7) 

Iai = J dV£, c V cP V a £ b V b £ a (A8) 

Ia 3 = J dVH c V c pli a V a (V-0- (A9) 

The integral Iai now has two derivatives acting on the pressure. Only one more integration by parts 
will be necessary to provide the term involving three covariant derivatives of pressure. The last integration 
by parts gives 

Iai = -J dVe^a (£TV c V b p) (A10) 

= -Jdv(v-z e^cVbP + rv„r c b v 6 v cP + re b rv a v b v cP ) . (ah) 

Notice that the second term is just the negative of the integral Iai- This allows us to rearrange the last 
equation to read 

Iai = -\ f dV ((V • £) £ 6 £ c V c V 6 p + C£ b rV a V h V cP ) . (A12) 

This step gives us the term involving three derivatives of the pressure. 

The integral Ib involves two derivatives of the displacement, which can be simplified through an inte- 
gration by parts to 



Ib = J dVV c £ b Va (p£°Vi£ c ) (A13) 
= J dV (V • (pOV„£ 6 V 6 £° + pVc^C V a V b r) • (A14) 



The last term in this expression is just the negative of 1b- Hence 

I B = ljdVV- (pO^a^bC- (A15) 

The final task is to transform Ic into a more useful form. This will be done by integrating by parts 
again, giving 

I c = J dVV ■ £V 6 ( P CV a e) = lex + I C 2 + I C3 , (A16) 

where 



lex = j dVV-i V bP Z a V a e (A17) 
I C 2 = J dVp V • a a V„(V • (A18) 
Ics = J dVp\7^V a eV b C. (A19) 

The final integration by parts allows us to rewrite Ici in the form 

lex = - j dVH b V a ((V • 0£°V 6 p) (A20) 

= --/^(e-vpe-v(v-e) + (v-o 2 ?-v P + (v-e)re b v a v b p). (A21) 

Putting together the results for I a, Ib and Ic, the original integral is 

I = - \ J dV (e^ c V a V fc \7 cP + 3(V • eCVrVbP - 3V ■ (p£) v a e b v 6 o 

- y ((v • c) 2 £ • Vp - p (v • o e a v a (v • o) (A22) 

While equation (A22) may not appear to be simpler than equation (Al), it will prove to be useful in 
simplifying the third-order Hamiltonian coupling r-modes. 

B. Order of Magnitudes for r-modes 

In the slow-rotation approximation, defined by Cl <C 1, the fluid displacement vector for r-modes has 
components which are of order 

C = 0(ti 2 ),t e = 0(l),^ = 0(l). (Bl) 

In the slow rotation approximation, the star's pressure deviates from spherical symmetry at second order in 
angular velocity, so that derivatives of the pressure are of order 

§- r P = Od),§- e P = O&),±p = 0. (B2) 

It then follows that 

TVaP = CdaP = 0(Q 2 ) (B3) 
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for r-modes. 

In equation (16) for the third-order r-mode Hamiltonian, a term involving two covariant derivatives of 
pressure appears. In terms of partial derivatives and Christoffel symbols for spherical coordinates, this term 
is 

re b v Q v fcP - ee(dad bP - r c ab d cP ). (B4) 

The non-vanishing Christoffel symbols for flat space written in spherical coordinates are 

r 2 Ke = r 2 T% = = = ^ T % = cot0 - ( B5 ) 

The leading order behaviour of equation (B4) is then 

e?V a \7 bP = U rP (£% + £%) = 0(1). (B6) 
Similarly, a longer but straight-forward calculation shows that 

V a £"V^ a = 0(1). (B7) 
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Table 1. Non- resonant Direct Couplings with an Unstable I = \m\ = 2 r-mode. The stellar model is a 
N = 1 polytrope with Ti = 1.9. In these tables, I and in refer to spherical harmonic angular quantum 
numbers and k refers the number of radial nodes in the eigenfunction. 



Mode A Mode B Mode C 



1 m k 1 m k 1 m k kabc Aojabc/Q 



First Generation Couplings 



2 


2 





2 


2 





5 


-4 





2.01c-02 


1.07 


2 


2 





2 


2 





5 


-4 


1 


1.54e-02 


1.07 


2 


2 





2 


2 





5 


-4 


2 


1.33e-02 


1.07 


2 


2 





2 


2 





5 


-4 


3 


1.20c-02 


1.07 


2 


2 





2 


2 





5 


-4 


4 


l.lle-02 


1.07 


2 


2 





2 


2 





5 


-4 


5 


1.04c-02 


1.07 


2 


2 





2 


2 





5 


-4 


6 


9.69e-03 


1.07 


2 


2 





2 


2 





5 


-4 


7 


9.58e-03 


1.07 


2 


2 





2 


2 





5 


-4 


8 


9.25e-03 


1.07 


2 


2 





2 


2 





5 


-4 


9 


9.06c-03 


1.07 


Second Generation Couplings 


2 


2 





5 


-4 





4 


2 





-8.10c-02 


0.60 


2 


2 





5 


-4 





6 


2 





2.33e-01 


0.49 


2 


2 





5 


-4 





8 


2 





9.83c-04 


0.45 


2 


2 





5 


4 





6 


-6 





-9.68e-02 


0.91 


2 


2 





5 


4 





8 


-6 





3.05c-02 


0.92 



Table 2. Comparison of coupling coefficients between two I — m = 2 r-modes and the I — 5, m = —4 

r-modc for different polytropic equations of state. 



N 


7=1 + 1/^ 


I T x 


i - 7/n 


K 


3 


4/3 


1.2 


-1/9 


-3.00e-03 


2 


3/2 


1 


-1/2 


-1.21e-02 




2 


1.99 


-5e-03 


2.16e-02 




2 


1.9 


-0.05 


2.01e-02 




2 


2.1 


0.05 


2.35e-02 




2 


2.01 


5e-03 


2.20e-02 




2 


2.0008 


4c-04 


4.00c-02 
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Table 3. Largest Coupling Coefficients. The stellar model is a N = 1 polytrope with T\ = 1.9. 



Mode A Mode B Mode C 



1 


m 


k 


1 


m 


k 


1 


m 


k 


KABC 




2 


2 





8 


-1 





9 


-1 





1.21c+01 


0.62 


2 


2 





8 


-1 


1 


9 


-1 


1 


1.07c+01 


0.62 


2 


2 





9 


-1 





10 


-1 





1.93e+01 


0.63 


2 


2 





9 


-1 


1 


10 


-1 


1 


1.71c+01 


0.63 


2 


2 





10 


-1 





11 


-1 





2.92c+01 


0.63 



Table 4. Indirect Resonant Couplings. The stellar model is a TV = 1 polytrope with Ti = 1.9. 



Mode A Mode B Mode C 



1 


m 


k 


1 


m 


k 


1 


m 


k 


KABC 


5 


-4 





5 


2 





5 


2 





-1.35e-02 


5 


-4 





5 


2 





5 


2 


1 


-1.46e-01 


5 


-2 


1 


5 


1 





5 


1 





1.40e+00 


7 


-2 





7 


1 





7 


1 


1 


-2.08c+00 


9 


-2 


1 


9 


1 





9 


1 





1.03e+01 
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Fig. 1. — Plot of the detuning versus m for values of I < 10. The detuning is defined in equation (40). 



